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1 Introduction 

Vector boson pair production is an outstandingly important process at high energy hadron 
colliders. Its measurement allows precision studies of the electroweak interaction, thereby 
testing in detail the SU {2)lxU{1)y gauge structure and the matter content of the Standard 
Model of particle physics. The various combinations of vector boson pairs [ZZ, W~'^W~, 
77 , ZW^, Z 7 , lead to spectacular final state signatures (leptons, photons, missing 

energy), that are often equally relevant to searches for new physics or studies of the Higgs 
boson. The Higgs boson decay into two vector bosons is among the cleanest signatures for 
Higgs production, and offers a broad spectrum of observables. 

Precision studies of the electroweak interaction often focus on the pair production of on- 
shell gauge bosons, while new physics searches and Higgs boson studies precisely veto these 
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on-shell contributions, such that the remaining background processes are dominated by off- 
shell gauge boson pair production. For both on-shell and off-shell production processes, it is 
therefore very important to have a precise prediction of the Standard Model contributions, 
in order to match the anticipated experimental accuracy of measurements at the LHC, 
which is usually in the per-cent range. At this level of precision, next-to-leading order 
(NLO) corrections in the electroweak theory and next-to-next-to-leading order (NNLO) 
corrections in QCD are indispensable. 

For all vector boson pair production processes, NLO QCD corrections [1-6] as well as 
large parts of the NLO electroweak corrections [7-14] are available. These calculations are 
fully differential in all kinematical variables, and usually include the leptonic decays of the 
vector bosons. The derivation of NNLO QCD corrections to vector boson pair production 
can build upon calculational techniques [15, 16] that were originally developed for the 
Drell-Yan process [17, 18] or for Higgs boson production in gluon fusion [15, 16], which 
have the same QCD structure due to their colour-neutral final state. As a new ingredient, 
each vector boson pair production process at NNLO requires the two-loop corrections to 
the basic scattering amplitude for quark-antiquark annihilation: —)• H 1 V 2 . These have 

been known for a while already for 77 [19, 20] and H 7 [21, 22] production, enabling the 
calculations of these processes [23, 24] to NNLO accuracy. 

Compared to the above, the two-loop matrix elements for the production of a pair 
of massive vector bosons require a new class of two-loop Feynman integrals: two-loop 
four-point functions with massless internal propagators and two massive external legs. 
Recently, very important progress has been made on these. For the case of equal vector 
boson mass, these integrals were derived in [25, 26], and used subsequently to compute the 
NNLO corrections to the on-shell production of ZZ [27] and W^W~ [28]. The integrals 
for the most general case of non-equal masses were derived in [29-31], which allowed to 
construct the full two-loop helicity amplitude for qq' —)• V 1 V 2 in [32]. A subset of these 
integrals was derived independently in [33, 34] and used in the derivation of the fermionic 
NNLO corrections to 7 * 7 * production [34]. In this paper, we perform an independent 
rederivation of these integrals and optimise our solutions for numerical performance. They 
are used subsequently for a validation of the two-loop helicity amplitudes of [32] , uncovering 
an error in their original results. We present a public implementation for the numerical 
evaluation of these amplitudes, which in the future will allow the calculation of NNLO 
QCD corrections to arbitrary electroweak four-fermion production processes. 

The paper is structured as follows: in Section 2, we introduce the partonic current 
for vector boson pair production and describe its decomposition into Lorentz structures. 
Taking into account the vector boson decays into leptons, we present the helicity amplitudes 
for four particle final state in Section 3. A detailed description of the calculation and the 
different contributions to the amplitude is given in Section 4. The computation of the 
master integrals and their optimisation is presented in 5. In Section 6 , we describe the 
subtraction of UV and IR counter terms, and in Section 7 we list the numerous checks 
we performed on our results. In Section 8 we present our C++ implementation for the 
numerical evaluation of the amplitudes and use it to produce numerical results. Finally, we 
conclude in Section 9. In Appendix A, we document the interference of the two-loop and 
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tree amplitudes for the production of on-shell vector boson pairs, which was used in the 
calculation of the NNLO corrections to pp —?■ ZZ [27] and pp —?• WW [28]. Appendix B 
contains the derivation of Schouten identities for the leptonic amplitudes, and Appendix C 
describes the conversion of our results between different schemes for the subtraction of 
infrared singularities. We provide computer readable files for our analytical results and 
our C++ code for the numerical evaluation of the amplitude on our VVamp project page on 
HepForge at http://vvamp.hepforge.org. 

2 Lorentz structure of the partonic current for qq' — )■ V 1 V 2 

Let us consider the production of two massive electroweak vector bosons in qq' annihilation: 

q{Pi) + 4 {P2) —^ Vi{p3) + V 2 {pa) (2.1) 

with 

P?=pi = 0, pI^O, p2^0, (2.2) 

where the two vector bosons are off-shell and V 1 V 2 = ZZ, W~^W~, 77 , ZW^, Z^, 

We dehne the usual Mandelstam variables 

'S = (P 1 +P 2 )^, ^=(Pl-P3)^ U=ip2-P3f, (2.3) 


such that 

s 1 u = P3 P4 . 

The physical region of phase space is bounded by tu = p|p| such that 

s > (+ Y^) > \{pI+pI- s- k) <t<]^{pl+pl-s +k) (2.4) 
where k is the Kallen function 

{s,pi,pl) = \/s‘^+pi+pi- 2 {spl +pIpI+pIs) . (2.5) 

Let us consider the partonic amplitude for the production of the two off-shell vector 
bosons L 1 V 2 

S{s,t,pl,pl) = S^^’'{pi,p2,P3)e^{p3T el{pA)* , ( 2 . 6 ) 

where 63 and £4 are the two polarisation vectors of Vi and V 2 respectively. In this notation, 
we keep an overall factor implicit, where e is the positron charge. 

In order to calculate the partonic current {pi,p 2 -,P 3 ), we consider its tensorial 
structure. Lorentz invariance restricts it to be a linear combination of 17 independent 
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structures 


S^‘'ipl,P2,P3) = u{p2) J^3u{pi) [FiKPi +F 2 P 1 P 2 +-^3PiP3] 

+ u{p2) ^ 3 uiPi) [ FaP2P1 + F3 P2P2 + FeP2P3 ] 

+ u{p2) ^ 3 u{pi) [ F7 p^Pi + Fs p'^P2 + Fg p'^p'^ ] 

+ u{P 2 ) l^uipi) [ -^10 Pi + FiiP2 + Fi 2 P3 ] 

+ u{p2) l''u{pi) [Fi3P^ + Fi 4P^ + Fisp^] 

+ u{p2) l^'^3l''u{pi) Fie 

+ u{p 2 )Y'^ 3 'y^u{pi) Fn, (2.7) 

where the form factors Fi ,..., Fn are scalar functions of the Mandelstam variables s, t,p\,p\ 
and of the number of space-time dimensions d. To further constrain we choose the 
Landau gauge for the electroweak vector bosons with the transversality condition 


es • P3 = £4 • P4 = 0 , 


and the sum over polarisations 

E( 4)'4 = -s"" + #. 

pol ^3 

^ ■ 

pol ^4 


( 2 . 8 ) 


(2.9) 


Imposing condition (2.8) we can reduce the number of independent tensor structures 
to ten [35, 36], which can be chosen as 


= u{p2) h^iPi) PiPi , 
= u{p2) hu{pi) P 2 P 1 , 
= HP 2 ) l^u{pi) p'( , 

= u{p2)-f‘'u{pi)p>^, 

Tg" = u{p2) l^'^3l''u{pi) , 


T!^'" = u{p2)hu{Pl)PiP2 , 
Fi'" = u{P2) hu{pi) P 2 P 2 > 
Tr= hp 2) i^uipi) P 2 , 

= HP2)Yu{pi)P2^ 
Tig = U{P2) Yhl^'^iPl) ■ 


( 2 . 10 ) 


Without any loss of generality we can thus write the partonic current as 

10 

S^''{pi,P2,P3) = Aj{s,t,pl,pl)Tlf'', ( 2 . 11 ) 

i=i 


where we introduced the new physical form factors ^di,..., j4io, which are again scalar 
functions of the Mandelstam variables s,t,pg,p‘l and of the dimension d. 

Note that in deriving (2.11) no assumption has been made on the dimensionality d, 
such that this decomposition is valid for any continuous values of d. Its structure has 
been constrained using solely Lorentz and gauge invariance and is therefore true at every 
order in perturbation theory. On the other hand, the scalar coefficients Aj{s,t,pg,p1) 


- 4 - 




contain the explicit dependence on the perturbative order at which they are computed. 
These coefficients can be extracted from the amplitude by applying appropriate projecting 
operators on the latter. The projectors themselves can be expanded in the same tensorial 
basis: 

10 

= E ’ i = 1,10 , (2.12) 

i=l 

where the coefficients Bji are functions of the Mandelstam variables and of the 

dimension d. They can be determined imposing that 


^ [es^, £4^. eLi eL. ] = A ,. (2.13) 

pol 

Note that the contraction is performed in d dimensions and at every stage one should 
always recall to use the polarisation sum in (2.9). For later convenience we introduce also 
the following scalar quantities: 

T = E [esMi U,. eLi 4.2 ] > (2-14) 

pol 

which are related to the coefficients Aj according to 

10 

= (2.15) 

i=l 


with the same coefficients Bji as in (2.12). These quantities (rather than the coefficients Aj) 
are particularly useful in order to build up the contractions of the n-loop amplitudes with 
the tree-level ones (see Appendix A). We provide explicit expressions for Bji in computer 
readable format on HepForge. 

The partonic current receives contributions from QCD radiative corrections and can 
be decomposed perturbatively as 

Siiu{pi,P2,P3) = Slj^J{pi,P2,P3) + SAJ{pi,P2,P3) + Sj^J{pi,P2,P3) + 0(0^) . 

(2.16) 

Obviously also the un-renormalised tensor coefficients Aj (or, equivalently the Tj) have the 
same perturbative expansion of the partonic amplitude 


II 

+ 

/as' 

) 

V27r/ 

1 3 

(0) 

p - Tj 

+ ( 

^a^N 

- 

K27r) 

3 


+ 


+ 


27r/ 
(£)^ 


4(2)+0(a3), 
+ C’(as), 


j 

( 2 ) 


where the dependence on the Mandelstam variables is again implicit. 


(2.17) 

(2.18) 


3 Helicity amplitudes for gg' — )■ I/ 1 V 2 —)■ 4 leptons 
In physical applications we are interested in the processes 

q{pi) + q {P 2 ) Vi{p^) + V' 2 (p 4 ) h{Pb) + kiPe) + hip?) + kiPs) (3.1) 


where each of the two off-shell electroweak vector bosons can decay to pairs of leptons, 
such that p 3 = P 5 + P6 and p^ = P 7 + Ps- Let us first focus on the general structure of the 
helicity amplitudes for this process. Schematically these amplitudes can be written as the 
product of the partonic current and two leptonic currents mediated by the 

propagators of the two off-shell vector bosons P^^{q) 

M(P5,P6,P7,P8;P1,P2) = S>^’'{pi,P2,P3) P^piP3)Lp{P5,P6) PuHP 4:) La{p7,P8) , (3.2) 


where we stripped off electroweak couplings not relevant here and postpone their discussion 
to the presentation of the full amplitude in (3.18) below. In the ii^-gauges the propagator 
of a vector boson V reads 



Dv{q) 


(3.3) 


with 


A 


iq,0 = + {1-0 


q^qu \ 

’ 


(3.4) 


D-y*{q) = q‘^, Dz,w{q) = {q^ - my+ iTvmv) , (3.5) 

where my is the mass of the gauge boson and Ly is its decay width. While the Landau 
gauge used in the previous Section corresponds to ^ 0, the term proportional to (1 — 0 

effectively vanishes for any ^ since the electroweak vector bosons are directly coupled to 
massless fermion lines. 

By fixing the helicities of the incoming partons and of the outgoing leptons one sees 
that the left- and right-handed partonic production currents can be written as 

S^L{pi,pt,P3) = v+{p2)'^^''u-{pi) = (2|r'^'"| 1] , (3.6) 

^^(P^Py^Ps) = u-(p2)r'^"n+(pi) = [2 1 ) , (3.7) 

where the are rank two-tensors and contain an odd number of y-matrices. We note in 
passing that, by complex conjugation, one gets 

[S>^^{pt,P 2 .pOV = ([2 1) )* = (2 1] = S>l0pi.pt,p0 , for all . 

The left- and right-handed leptonic decay currents, on the other hand, can be written as 

pLiPb^Pe) = u-{Pb)l^v+{pO = [6|7^|5) = (5|7^|6], (3.8) 

pRipt^Pb) = u+{p0i^V-{p0 = [Sly'^lG) = (L^(P5 ,P6 ))* = ^^(Pe ,p^). (3.9) 
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Note in particular that, as far as the lepton currents are concerned, a permutation of 
the external momenta corresponds to a flip of the helicity. All possible helicity amplitudes 
can be therefore obtained from the two basic amplitudes 

MlLl(P1,P2;P5,P6,P7,P8) = {pT :Pt ^Ps) LLfj,{P5 ,pt)LLu{P7 ,Ps) , (3-10) 

^RLl{P1,P2;P5,P6,P7,P8) = {pf ,P2 ,P3) ,pt)LLu{P7 ,Ps) , (3-11) 

by simple permutations of the leptonic momenta. In particular we find 

Mllr{P1,P2]P5,P6,P7,P8) = M.lll{P1,P2;P5,P6,P8,P7) , 
'!^LRl{P1,P2]P5,P6,P7,P8) = ^LLl{P1,P2;P6,P5,P7,P8) , 
'^LRr{P1,P2]P5,P6,P7,P8) = Mlll{P1, P2] P6,P5, P8, P 7 ) , 
'^RLr{P1,P2]P5,P6,P7,P8) = ^RLl{P 1, P2] Pb, P6, P8, P 7 ) , 
'^RRl{P1,P2]P5,P6,P7,P8) = Mrll{P 1, P2-, P6, Pb, P7, Ps) , 
'^RRR{Pl,P2;Pb,P6,P7,P8) = M.RLL{pi,P2;P6,Pb,P8,P7) ■ (3.12) 

In order to put together the helicity amplitudes in their final form we need also to take 
into account the electroweak couplings of the gauge bosons to the partonic- and leptonic- 
currents which we have been kept implicit so far. We follow [37] and parametrise the 
coupling of a vector boson 1/ to a fermion pair / 1/2 as 


= z e, where e = V^ira is the positron charge , (3.13) 

such that all fermion charges are expressed in units of e and 

W* = ■ (3.14) 

The left- and right-handed interactions are equal for a purely vectorial interaction. De¬ 
pending on the different kinds of gauge bosons we have 


L 


/1/2 “ 


z _ I 3 - sin e^ef, ^ 

sin0,„cos0,„ ’ 


^flf2 1 

z single/, ^ 

cos9,„ ’ 


tW _ 


1 




</ 2 = 0 , 


(3.15) 

(3.16) 

(3.17) 


where again the charges are measured in terms of the fundamental electric charge e > 0 
and e/j /2 is unity for /i / / 2 , but belonging to the same isospin doublet and respecting 
charge conservation, and zero otherwise. 

Putting everything together we find for the two independent helicity amplitudes for 
qq' —> V1V2 —> hhhh 


-^ALL^(P1’P2;P5,P6,P7,P8) = (4770;)^ 


tVi r V2 


DviiP3)Dv2{P4) 


MxLL{Pl,P2]Pb,P6,P7,P8) , (3.18) 
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where X = L,R and we have bracketed out the tree-level dependence on the electric charge 
i(47rQ;)^ and on the couplings with the decay lepton currents. Obviously the correspond¬ 
ing helicity amplitudes for right-handed leptonic currents can be obtained by the simple 
exchange -H- R^-j, together with pi -H- pj. 

Once the tensor structure (2.10) is given, we can perform the contraction with the 
leptonic decay currents and fix the helicities of the incoming and outgoing fermions. This 
enables us to cast the two independent helicity amplitudes Mm and Mrll in the familiar 
spinor-helicity notation [4, 38]. In doing so, one assumes that the external states are 4- 
dimensional and this allows to obtain one further Schouten identity between the 10 tensors 
structures, such that one ends up with 9 independent form factors. Our derivation is spelled 
out in detail in Appendix B. As a result, we obtain 

Mlll{pi,P2;P5,P6,P7,P8) = [12^3 2) ^Ei (15)(17)[16][18] 

+ E 2 (15) (27) [16] [28] + Es (25) (17) [26] [18] 

+ Ei (25) (27) [26] [28] + E^ (57) [68]} 

+ Eg (15) (27) [16] [18] + Ej (25) (27) [26] [18] 

+ Es (25) (17) [16] [18] + Eg (25) (27) [16] [28] , (3.19) 

MiiLLiPi,P2;P5,P6,P7,P8) = [22^3 1) (15)(17)[16][18] 

+ E 2 (15) (27) [16] [28] + E 3 (25) (17) [26] [18] 

+ ^4 (25)(27)[26][28] + .E5(57)[68]} 

+ Eg (15) (17) [16] [28] + Ej (25) (17) [26] [28] 

+ Eg (15) (17) [26] [18] + Eg (15) (27) [26] [28] , (3.20) 

where 

[1^32) = [15](52) + [16](62), [ 2 ^ 3 !) = [25](51) + [26](61) , 

and the 9 form factors Ej are linear combinations of the form factors Aj 


El — Ai, 

s 

2 

E2 = A2 + - {Ag — Aio) , 
s 

s 

2 

E3 = Ag -(Ag - Alo) , 

s 

2 

Eg = 2 Ag — - [(u — S — P3)Ag + {t — P4 )Aio] , 

II 

2 

Eg = 2 Ag — - [(t — S — p|)Aio + {u — P4)A9] . 

Eg — 2 (Ag -I- Alo) . 

(3.21) 


In the following, we will consider a perturbative expansion of the form factors Ej defined 
completely analogous to that of the coefficients Aj in (2.17). We note that the expres¬ 
sions (3.19) and (3.20) are formally identical to the corresponding formulas derived in [32], 
such that our form factors Ej can be mapped one to one to the Fj defined in [32]. 
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Figure 1. Representative Feynman diagrams for classes B, C and Fy relevant for the produc¬ 
tion of two electroweak vector bosons at the two-loop level. All of these classes receive contributions 
both from planar and non-planar diagrams. 


In the next Section we will describe how to compute form factors Ai,...,Aio and 
therefore also form factors Ei,... at tree-level, one-loop and two-loop order, follow¬ 
ing a straightforward diagrammatic approach. In particular, we will discuss the different 
electroweak coupling arrangements C contributing to the functions Aj and Ej, 


A.-F ■ 

Q 

j = l,...,10, 

-p _ X \ ^ p)^jAV 2,[C] pic] 

Ej — E- , 

j = 1,... ,9, 


where denotes a coupling factor, A is the helicity of the incoming quark, and ii, 

Z 2 are the colours of the incoming quark and anti-quark, respectively. 

We want to stress once more an important point. Reducing the 10 coefficients Aj to 
the 9 coefficients Ej required the assumption that the external states can be treated as 
4-dimensional. In order to avoid any loss of information, we will work considering the Aj 
as fundamental objects (derived in d dimensions throughout) and refer to formulas (3.21) 
in order to reconstruct the Ej explicitly. 

4 Organisation of the calculation 

The calculation of the two-loop helicity amplitudes can be set up in a way that is indepen¬ 
dent on the nature of the vector bosons considered, by organising the Feynman diagrams 
contributing to any such process into different classes. We find in particular that, as long as 
we limit ourselves to QCD corrections only, at any given number of loops, seven different 
types of diagrams can contribute, depending on the arrangement of the external vector 
bosons. 

Class A collects all those diagrams where both vector bosons are attached on the external 
fermion line, such that Vi is adjacent to the quark q{pi). In the case of a left- 


- 9 - 















handed (right-handed) quark amplitude these diagrams are proportional to 

( 7?') 
y^qq" ^q"q')- 


Class B collects all those diagrams where both vector bosons are attached on the external 
fermion line, such that Vi is adjacent to the antiquark q'{p 2 )- Also these diagrams, 
in the case of a left-handed (right-handed) quark amplitude, are proportional to 


J-Vl r V2 /Dr] 
^q'q" ^q^'q '^^q'q' 


Vl nV2 




Class C contains instead all diagrams where both vector bosons are attached to a fermion 
loop. These diagrams are proportional to the charge weighted sum of the quark 
flavours, which we denote as A^ViVi; depending on nature of the final state bosons. 
In the general case, these diagrams yield two different contributions. In the first one, 
which is proportional to the sum of the vector-vector and the axial-axial couplings, 
all dependence on 75 cancels out. The vector-axial contribution, instead, is linear 
in 75 . Nevertheless, this last contribution is expected to always vanish identically 
for massless quarks running in the loops, for any choice of Vi and V 2 , due to charge 
parity conservation [32, 39-41]. Taking this into account we find 


"77 = IE[K*)'+W«.)' 

i 


"^7 = 5 E > 

i 

^ww = ’ (^- 1 ) 


where the indices i, j run over the flavours of the quarks in the loop and Lq-q- = Rq-q- 
Of course, ^ • e^. and, due to charge conservation, = Nwz = 0. 


Class Di contains all diagrams where Vi is attached to a fermion loop and V 2 to the 
external fermion line. Up to two loops, the diagrams in this class must sum up to 
zero due to Furry’s theorem. 


Class D 2 contains all diagrams where V 2 is attached to a fermion loop and Vi to the 
external fermion line. At two loops the diagrams in this class, as for the previous 
case, must sum up to zero due to Furry’s theorem. 


Class E contains all diagrams there Vi and V 2 are attached to two different fermion loops. 
These diagrams contribute only starting from three-loop order and we can ignore 
them. 


Classes Fv collect the form-factor diagrams where the production of the two vector 
bosons Vi,V 2 is mediated by the exchange of another vector boson V. Depend¬ 
ing on the type of vector bosons Vi , V 2 there can be more than one such class due 
to different intermediate vector bosons. In the case of a left-handed (right-handed) 
quark amplitude these diagrams are proportional to L'^^, CVV 1 V 2 cy V 1 V 2 )) where 
CVV 1 V 2 is the electroweak coupling of the triple gauge boson vertex defined for all 
particles and momenta outgoing as 

'^v7v 2 b,c) = iecv V 1 V 2 [ (« - bfg^P + {b - cYg^^'' + {c- aYg"'^ ] (4.2) 
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with 


— Cvi/=F-yVy± — C\Y±y/^^ — ±1 , 

cziyiVFT = cvj/=fzvf± = — Tcot0^ . 


(4.3) 


It is clear that, depending on the nature of the vector bosons Vi, V 2 and on the loop 
order, not all classes above will give non-zero contribution. At tree-level, for example, only 
classes A, B and Fy can contribute. The same is true also at one loop, provided that 
we limit ourselves to QCD corrections only. The situation changes at two loops, where 
also diagrams for classes C, Di and D 2 occur. Notice moreover that the form-factor type 
diagrams in class Fy are relevant only for the production of IT 7 , WZ or WW pairs. 

Up to two loops, we can thus restrict the summation in (3.22) to C = A, B, C, Fy. We 
show representative diagrams in Figure 1. For the coupling factors we have 


qL,ViV2,[A] 

^qq' 

qL,ViV2,IB] 

^qq' 

QLyiV2,[C] 

^qq' 

qL,ViV2,[Fv] 

'^qq' 


tVi tV2 

q q// ? 

^RyiV2,[A] 

^qq' 

tVi tV2 

^qfqff ^q"q ’ 

^RyiV2m 

^qq' 

Ny^y^Sqq! , 

QRyiV2,[C] 

^qq' 

Lqq>Cvv^V2 

s — rriy — iTymy ’ 

QRyiV2\Fv] 

^qq' 


— ^^q"q' ’ 

_ pVl pV 2 

— , 

= Ny^y^dqqi , 

_ Rqq'CVV^V2 
s — rriy — iVymy 


(4.4) 


With these dehnitions, the value of the coefficients do not depend on the nature 

of the mediating vector boson V, such that in particular 


j^F^],{n) _ j^Fz\,{n) _ ^Tw],(«) 




(4.5) 


We have computed the coefficients Aj for the different classes of diagrams contribut¬ 
ing at tree level, one loop and two loops, namely with C = 

A, B, C, Di,D2, F. 

At tree-level order we find 


A^A,(o) 

2 

.[A,(o) 

^10 

1 

’ 

4[A,(0) 

= 0 , 

j = 1 

, ..., 6 , 

8,9, 

aFUo) 

2 

— +“ ) 
u 

4[i?],(0) 

1 

? 

u 


= 0 , 

j = 1 

,...,7, 

10 , 

aT 1 -(o) 

= = +2 , 

4[C,(0) 

^9 

_ 4^140) _ -1 

— — i , 

^[C,(o) 

= 0 , 

j = 1 

,...,6 

• (4.6) 


We can notice immediately that, as far as the form-factor type diagrams are concerned, 
any n-loop QCD corrections will not modify the structure of (4.6), and in particular we 
have 




[C,( 0 ) 


where are the n-loop QCD corrections to the quark form-factor. 


(4.7) 
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Let us discuss the features of our Ej set of coefficients, which is relevant for the four¬ 
dimensional helicity amplitudes for the full 2—7-4 process. We consider crossings of external 
legs described by the permutations 


-^12 ■■= Pi P2 ^ {i u] 

vr34 := P3 P4 => { t ^ M O P4 } . 

4C] 


(4.8) 


and focus on the behaviour of the for the non-trivial cases C = A, B, C. From the 
exchange of quark and anti-quark, 7ri2 we find for the amplitudes 


“ill = « P 2 ), mSi = -Mgi„(pi « n ), 

from which one can directly obtain 
E^f^{s,t,pl,pl) = -E'f\s,u,pl,pl), 


E'f^{s,t,pl.,pl) = -Ef\s,u,pl,pl), 
= -Ef\s,u,pl,pl), 
E'f^isA^pl^pl) = -E'f\s,u,pl,pl), 
E'f\s,t,pl:Pl) = -Ef\s,u,pl,pl), 
E'f\s,t,pl,pl) = -E\j^\s,u,pI,pI) , 
E'f\s,t,pl,pl) = -Ef\s,u,pl,pl), 


Ef^{s,t,P%pl) = -Ef\s,u,pl,pl), 
E'f^isA^phpl) = -E'f\s,u,pl,pl), 
E'PisA^phpl) = -Ef\s,u,pl,pl), 
Ef^{s,t,pl,pl) = -Ef\s,u,pl,pl), 
Ef^{s,t,pl,pl) = -Ef\s,u,pl,pl), 
Ef^{s,t,pl,pl) = -Ef\s,u,pl,pl), 
Ef^{s,t,pl,pl) = -Ef\s,u,pl,pl) , 


From exchange of the external vector bosons, 7r34, we have 


M 

^XLL = ^^xLiP^ ^ Pi) : 

) 

mSl = 

Mf4(P3 

^P4), 

with A = 

L,R, 

which implies 









is,t,pl,pl) 

= -EPis, 

u. 

^pIps), 

4'^4: 

Lpi,P4) 

= +4''4, 

u,pI,P3) . 


{s,t,pl,pl) 

= -4^4, 

u. 

>pI,P3)^ 

Ef\s., 

Lpi,P4) 

= -Efhs, 

w,pipi) , 


{s,t,pl,pl) 


u. 

,pI,P3) , 


Lpi,pi 

= -eP{s, 

^^,pipi) , 


{s,t,pl,pl) 

= -E^is, 

u. 

<pI,P3) , 

E^'Hs., 

Lpi,pi 

= -Ef^\s, 

^^,pipi) , 

4"’ 

{s,t,pl,pl) 

= -4'"(^, 

u. 

<pI,P3) 2 

eP{s.. 

Lpi,pi 

= -EPis, 

^i,pipi) , 


{s,t,pl,pl) 

= +4'"(^, 

u. 

7 pi 4) 2 

eP{S: 

Lpi,pi 

= +eP{s, 

^i,pipi) , 


{s,t,pl,pl) 

= +4'"(^2 

u. 

.pi pi) 2 

eP{s.^ 

Lpi,pi 

= +eP{s, 

^^,pipi) , 

4"’ 

{s,t,pl,pl) 

= +4'"(^, 

u. 

.pi pi) 2 






(4.9) 


(4.10) 


(4.11) 


(4.12) 

Similar but slightly more involved relations can be derived for the primary set of coef¬ 
ficients, Aj, but we don’t list them here for brevity. We have explicitly verified that 
relations (4.10),(4.12) for the Ej and the corresponding relations for the Aj hold for our 
results at tree level, one loop and two loops. 


While most of the coefficients Aj are zero at tree level, fewer of the Ej have this 
property. We find for class A 


^[A],(0) ^ 

0, 

eAUo) 

11;^ 

1 

eAua 


0, 

eAUA 

2 

^[A],(0) 

_ 

ni^-Ps) 
st ’ 

pAm 

n^t-pl) 
st ’ 


for class B 






0, 

Ewm 

2 

su ’ 

eIbUo) 


0, 

eAW) 

2 

u 

eAUA 


^{s-u + pD 

p[Bm 

^8 

^{s-u + pD 
su ’ 


su ’ 

and for class F 





0, 

eAW) 

= 0, 



0, 

E^P^A) 

= -4, 

pAm 


+4, 

E^mA 

= -4, 

pAUo) 


2 

st ’ 

^ {s-t + pl) 

st 


(4.13) 


2 

su' 

^ {u-pD 

su 

o (^-P4) 
SU ’ 


(4.14) 


0 , 

+4, 

-4. (4.15) 


As discussed above, class C contributions enter only at the two-loop level. 

The calculation of the coefficients Aj and thus Ej proceeds as follows. The diagrams 
belonging to class Ey are known [42]. They do not have to be recomputed and we will 
not refer to them anymore here. As far as the other classes are concerned, we produced 
all the tree-level, one-loop and two-loop Feynman diagrams with Qgraf [43]. The scalar 
coefficients Aj are evaluated analytically diagram by diagram by applying the projectors 
defined in (2.12) and summing over the polarisations of the external vector bosons as 
in (2.9). For the gluons we employ the Feynman-’tHooft gauge. All these manipulations 
are consistently performed in d dimensions. Upon doing this we obtain the coefficients in 
terms of a large number of scalar two-loop Feynman integrals. The latter are classified into 
three integral families, two planar and one non-planar. We have made use of Reduze 2 [44- 
47] in order to map all integrals to these integral families and to perform a full integration- 
by-parts reduction [48-51] of the latter to a set of master integrals. All intermediate 
algebraic manipulations on the Feynman diagrams have been performed using Form [52]. 
Once the coefficients Aj for the different classes of diagrams are known at the different loop 
orders, one can calculate the form factors Ej using (3.21). Since the expressions for the 
coefficients Aj (and equivalently those for the Ej) at two loops are very lengthy we decided 
not to include them explicitly in the text. Analytical expressions for the Aj, prior to UV 
renormalisation and IR subtraction, expressed as linear combinations of masters integrals 
and retaining full dependence on the dimensions d are available on our project page at 
HepForge. 
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5 Master integrals 

5.1 Computation via differential equations 


We computed all two-loop master integrals needed for our process with the method of 
differential equations [50, 53-55] and optimised the solutions for fast and precise numerical 
evaluations [26, 56, 57]. The master integrals for the case = p\ have first been calculated 
in [25, 26]. Here, we consider the case p\ ^ p\, for which the master integrals have been 
computed in [29, 30] for the first time. Our calculation provides an independent check of 
these results and improves them for numerical applications. In this Section we present our 
calculation and discuss qualitative aspects of the results. We provide the explicit solutions 
in computer readable format on HepForge. 

We find that all master integrals are described by the integral families presented in [26] 
for the case ^ p\ and crossings thereof. We start by determining a set of linearly inde¬ 
pendent master integrals for all relevant topologies using Reduze 2 [44]. For convenience, 
we stick to the normal form definitions for the master integrals given in [29, 30]. We 
supplement these definitions by new normal form definitions for eight factorisable topolo¬ 
gies corresponding to products of one-loop integrals. All our definitions are supplied in 
computer readable form on HepForge. 

We consider the master integrals of all integral families at the same time and eliminate 
multiple variants of equivalent master integrals using the shift-finder of Reduze 2. For 
this purpose we also identify crossed topologies and work out relations between crossed 
and uncrossed master integrals. Ignoring crossed variants and counting product topologies 
as two-loop topologies we find a total number of 84 independent master integrals. To 
apply the method of differential equations, we include also crossed versions for a couple of 
integrals, which appear in sub-topologies of non-planar topologies. In this way we assemble 
a minimal set of 111 master integrals suitable for the construction of a system of differential 
equations. 

We compute the partial derivatives of the master integrals with respect to all indepen¬ 
dent external invariants s, t, p^, p\ in terms of master integrals with the help of Reduze. 
The coefficients contain rational functions of the invariants and the Kallen function k, (2.5), 
associated to the two-body phase space. 

To rationalise the root k, we employ the parametrisation 


s = {1 + x)'^, p\ = — y^), 

t =-fh^x{{l + y){l + xy)-2zy{l + x)), = m^(l - , (5.1) 


(see eq. (2.9) of [30]). In this parametrisation, we define the vector of master 
M = (Mj), i = l,...,lll, using the integral measure 

2,2. f d^k dH 

Vl67r2 J ) y {2'nY (27r)'^ 

which absorbs the overall mass dimension m. Here, d = 4 — 2e and 


integrals 


(5.2) 


C'. 


(47r)" 


r(l + e) r 2 (l-e) 
F(1 - 2e) 


(5.3) 
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In the following, we will work directly in the physical region of phase space. Due to the 
specific choice of the master integrals [58, 59], the partial differential equations combine 
into the simple total differential. 


20 


dM (e; x,y,z) = A^d ln(4) M (e; x, y, z) 


(5.4) 


k=l 


where the matrices Ak contain just rational numbers and the alphabet is 


{Zi,..., / 20 } = {2, X, 1 + X, 1 - y, y, 1 + y, 1 - xy, 1 + xy, I- z,z, 


I + y - 2yz, 1 - y + 2yz, I + xy - 2xyz, 1 - xy + 2xyz, 

_ _ _9 _ _ _ _ _9 _ _ 

1 + y + xy + xy — 2yz — 2xyz, 1 + y — xy — xy — 2yz + 2xyz, 

_ _ _9 _ _ _ _ _9 _ _ 

1 — y — xy + xy + 2yz + 2xyz, 1 — y + xy — xy + 2yz — 2xyz, 

1 — 2y — xy + y'^ + 2xy‘^ — xy^ + dyz + 2xyz + 2xy^z, 

1 — y — 2xy + 2xy‘^ + — x^y^ + 2yz + 4xyz + 2x‘^y^z} . (5.5) 


Anticipating the solution, we included the letter 2 already, which is of course arbitrary 
at the level of the differential equations. While we found that it is possible to reduce 
the number of letters by forming appropriate ratios, a reduction of the alphabet is best 
performed using a different parametrisation, as we will see below. 

After expansion in e it is straight-forward to integrate the differential equations in 
terms of multiple polylogarithms 



(5.6) 


with G(0,..., 0; z) = ^ ln"'(z) for n zero weights and G(; z) = 1. For each order in e, we 
integrate the partial derivatives in z. This gives the solution up to a function of x and 
y. We employ the partial derivatives in x to determine this function, this time up to a 
function of y. Subsequent usage of the derivative in y fixes the boundary terms up to one 
constant per master integral for the given order in e. Despite the presence of nonlinearities 
in (5.5), the specific order of our integrations ensures that in fact only linear denominators 
occur in the respective integration variable. We integrate the master integrals through 
to weight 4, which corresponds to terms in the chosen normalisation. The necessary 
argument-change transformations for the multiple polylogarithms were derived using an 
in-house package, which employs fitting of constants using high precision samples obtained 
with [60]. 

In order to fix the integration constants, we consider the equal mass limit —)• 
which implies x —?• 1. This limit is smooth and our master integrals become simple linear 
combinations of the normal form integrals defined in [26], where the coefficients in this 
map are just rational numbers. We compute the limit at the level of our solutions and 
equate them to the real-valued solutions of [26]. Using the coproduct-augmented symbol 
calculus [56, 61-65], we find perfect agreement for all non-constant terms and easily fix the 
boundary constants of the present integrals. We also compared our results to the solutions 
of [29, 30] and hnd perfect agreement at the analytical level. 
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The solutions we obtained in this way are not ideal for our purposes yet, since their 
numerical evaluation is rather slow. Moreover, they contain spurious structures: the indi¬ 
vidual multiple polylogarithms contribute letters {1 — x, 1 -|- 2 + x + xy^, 1 -|- 2x -|-x^y^} 
which cancel for the master integral itself. In particular, the equal-virtuality limit x —?• 1 
is completely smooth as can be seen from (5.5), but the representation does not allow for 
an evaluation exactly in the equal-virtuality point. 

5.2 Optimisation of the functional basis 

We wish to cast our solutions to a new representation which allows for fast and stable 
numerical evaluations and is free of spurious letters. In order to achieve that goal we select 
a new basis of multiple polylogarithms where we do not force individual variables into 
the argument of G functions anymore. As a side effect, this gives us more freedom for a 
rational parametrisation, since we avoid problems due to non-linear denominators in the 
integration variable. It is convenient to choose new variables x, y, z and according to 

s = m^(l-|-x)(l-f xy), t = —rn?xz, p‘^ = m^, p\ = m?x‘^y (5.7) 

(see eq. (2.7) of [30]), which again rationalises the root n. We select the branch for which 
in the physical domain 

X > 0, 0 < y < z < 1, > 0. (5.8) 

This reparametrisation is actually not crucial for what will follow, but it decreases the num¬ 
ber of irreducible polynomial letters which will be convenient for our mapping procedure. 
Under crossings of external legs the parameters transform as 

7 ri2 : z ^ I + y - z (5.9) 

7r34 : z^l + y — z, x—)>l/(xy), rri^ ^ m^x‘^y. (5.10) 

In this parametrisation we factor out a normalisation of the form (5.2) but with m replaced 
by m. We find the alphabet 

{h,..., In} = {x,l + x,y,l - y,z,l - z,-y + z,l + y - z,l + xy,l + xz, xy -h z, 

o 

l-|-y-|-xy — z, l-|-x-|-xy — xz, 1 -|- y -|- 2xy — z + x yz, 

9 9 9 9 9 

2xy + X y + X y + z — x yz, 1 + x + y + xy + xy — z — xz — xyz, 

1 + y + xy + y"^ + xy"^ - z - yz - xyz} (5-11) 

at the level of the differential equations and also of the solutions through to weight 4. This 
alphabet is shorter than the previous one and can not be reduced further by forming ratios. 

We construct a new functional basis consisting of Li 2,2 functions, classical polyloga¬ 
rithms Li„ (n = 2, 3,4) and logarithms, similar to the approach taken in [26, 57]. The Li 2,2 
function can be written in G-function notation according to 

Li 2 , 2 (xi, X 2 ) = G fo, —, 0, —; . (5.12) 

V 3 : 1 x 2 J 
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Following the algorithm of [64] we generate functional arguments which are rational func¬ 
tions of X, y, z and do not lead to new spurious letters. This implies that the arguments 
factorise into the letters of our alphabet and their inverses. We can therefore systematically 
scan for admissable Li„ arguments by constructing power products of letters, their inverses 
and —1. A candidate argument xi is admissable exactly if 1 — xi factorises into the letters 
of our alphabet and —1, since only in that case the introduction of new letters is avoided. 
Admissable arguments for Li 2,2 functions are determined by forming pairs of admissable 
Li„ arguments and requiring for any such pair (xi, X 2 ) that the difference xi — X 2 factorises 
into the original letters and —1. 

For the amplitude we need to evaluate also independent master integrals with crossed 
kinematics, and we chose to implement these expressions explicitly for evaluation time 
optimisation purposes. We therefore directly construct a shared set of basis functions for 
uncrossed and crossed variants of the master integrals and consequently close our alpha¬ 
bet (5.11) under 7ri2 and by adding the letters 

{^ 18 , ^ 19 } = {-xy + z + xz + xyz, -y + z + yz + xyz}. (5.13) 


We require all functions to be single valued and real over the entire physical region of 
phase space. As in [26], we further tighten this constraint and select only those Ij\ 2 , 2 {xi,X 2 ) 
functions, for which their power series representation 


00 00 

Li2,2(2:1,X2) = X] X] 
11=112=1 


Ji 


( xiX 2)-^2 


{jl +j2f 


3 l 


(5.14) 


is convergent, that is, their arguments fulfil 


Xl| < 1, IX 1 X 2 I < 1. 


(5.15) 


We wish to express our master integrals in terms of these new functions and employ 
the coproduct-augmented symbol calculus for that mapping, see [56, 64, 65]. This step 
is computationally demanding due to the large number of possible candidate functions. 
Here, we profit from the reduction of the number of letters described above which leads 
to a smaller set of candidate integrals for a given maximal total degree of the arguments. 
Furthermore, we employ a particularly efficient technique for the symbol calculus, where 
we identify and match individual factors of products directly at the level of the symbol [66] . 
In particular, this means we never need to construct products of polylogarithms for our 
candidate functions which avoids a severe combinatorial blowup for the linear algebra 
routines. Using the coproduct we were able to express all master integrals in terms of 
our new set of functions described above. We stress that the success of this matching is 
not a priori obvious. The explicit solutions for all of the master integrals are provided on 
HepForge. 

Concerning our primary motivation for changing the functional basis, we observe that 
the new representation indeed allows for significantly faster numerical evaluations. For the 
numerical evaluation of the multiple polylogarithms we employ the implementation [60] 
in the GiNaC library [46]. The exact evaluation time and the speedup due to the new 
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functional basis depend on the chosen point in phase space and on the required precision 
of the result. We tested some samples and observed speedup-factors between 8 and 85 
when evaluating the 111 master integrals in the system of differential equations. For the 
benchmark point of [32], the numerical evaluation with default precision takes 2250 ms for 
the “traditional” G-functions (Section 5.1) and 120 ms for the “optimised” functions (this 
Section) on a single core of a standard desktop computer. 


6 UV renormalisation and IR subtraction 


Let us go back to the calculation of the helicity amplitude coefficients Aj (or equivalently 
the Ej). In order to simplify the notation for what follows we pick one of the form factors: 

n = Aj (or Ej) , for some j = 1,..., 10 (9), (6.1) 


in order to suppress the index j. The following discussions applies to any of the chosen 
form factors in the same way. 

We perform renormalisation of the UV divergences in the standard MS scheme which, 
in massless QCD, amounts to simply replacing the bare coupling ao with the renormalised 
one as = as(//^), where is the renormalisation scale. Since in our case the tree-level 
amplitudes do not contain any power of ag we require only the one-loop relation for the 
coupling 


'5*6 — ag /J 


1 


/3o 



( 6 . 2 ) 


where 


S'e = (47 r)^e , with the Euler-Mascheroni constant 7 = 0.5772..., (6.3) 


//q is the mass-parameter introduced in dimensional regularisation to maintain a dimen¬ 
sionless coupling in the bare QCD Lagrangian density, and finally /3o is the first order of 
the QCD /3-function 


/3o 


IIC*^ — AEp Nj 


with Ca = N , Cp 


-1 _ 1 
-, Ep — — . 

2iV ’ 2 


(6.4) 


We perform UV renormalisation at the scale = s, the invariant mass of the vector boson 
pair. Values of the helicity coefficients at different renormalisation scales can be recovered 
by using the renormalisation group equation. Since at a given loop order n the form 
factors are defined with all powers of the strong coupling factored out, the renormalised 
form factors are expressed in terms of the un-renormalised ones according to 


q{ 0) ^ Q(0),un ^ 

= S-^ , 

0(2) = 5,-2 57(2)’"“ _ ^ ^-1 J^(l),un _ 

After performing UV renormalisation, the amplitude contains residual IR singularities 
which will be cancelled analytically by those occurring in radiative processes at the same 


- 18 - 






order. Catani was the first to show how to organise the IR-pole structure up to two-loop in 
QCD [67]. In subtracting the poles from the one- and two-loop amplitudes we will follow 
a slightly modified scheme described in [68], which is better suited for the (^T-subtraction 
formalism. The two schemes are of course equivalent and we provide formulae to convert 
the results between the two schemes in Appendix C. We define the IR finite amplitudes at 
renormalisation scale in terms of the UV renormalised ones as follows 

5 ^( 2 ),finite ^ ^{ 2 ) _ j^(l) _ ^{ 0 ) ^ 


with 


hie) = irf\e) + I^ie), 


(6.7) 


/f^'(e) = - 




r(l-e) V s 


(4 + -+«Mc’p, = 


2\ 


1 in 


IT 


(3o 


hie) = -^hie)’^ + y [/i(2e) - /i(e)] + Kll<^f\2e) + H^e ), 

1 / 2\ 2e / (1) \ 

= 5 \n) +CF<ii. 

and the constants are defined as 
«' = 0 . 

qii 10^ ^ f 1214 67^ /164 5^\ , 

4’ - jC3/3o+ (-^ + ^C 2 j + (-^ - 42 ) Nf, 

/ 17 88 \ / 2 1 fi 

= (—3 -I- 24^2 — 48 C 3 ) Cp + ( -^C2 + 24^3 j Cp Ca+ ( 3 + “^^2 



( 6 . 8 ) 

(6.9) 

( 6 . 10 ) 


Cp Np. 
( 6 . 11 ) 


Note that in these equations all imaginary parts are already explicit prior to expansion in e. 
Setting fjp = s, we calculated the finite remainder of the Aj for e —?■ 0 in the (^T-subtraction 
scheme. We provide the explicit analytical results on our project page at HepForge. It is 
straight-forward to convert our finite results obtained in the gT-scheme to Catani’s original 
scheme, see Appendix C. 


7 Checks on the amplitudes 

We performed different checks on our amplitude, which we enumerate here. 
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1. First of all, we started off by computing the 10 form factors Aj of (2.11) for the 
different classes of diagrams C = A, B,C, Di, D 2 , and we explicitly verified that, 
according to Furry’s theorem, the diagrams in classes Di and D 2 independently sum 
up to zero. 

2. From the Aj we computed the 9 form factors Ej in (3.19) and (3.20), and we verihed 
that, both prior to as well as after subtraction of UV and IR poles, all symmetry rela¬ 
tions described in (4. 10), (4. 12) and the corresponding ones for the Aj, are identically 
satisfied. 

3. We verified that the poles of the one-loop and two-loop amplitudes are correctly 
reproduced by Catani’s formula [67], which provides a strong check on the calculation. 

4. For the NNLO computation of on-shell ZZ and W^W~ production [27, 28] we per¬ 
formed a dedicated calculation, directly for the squared amplitude, employing our 
equal-mass master integrals [25]. The tree and one-loop contributions have been 
found to agree with the analytical results of [69, 70] and with numerical samples 
obtained with OpenLoops [71]. Starting from our general results for the amplitude 
in the off-shell case, we re-derived the squared amplitudes for on-shell ZZ and WW 
production as described in Appendix A and found full agreement through to two- 
loops. 

5. We performed a thorough comparison of our results with an earlier calculation of 
the two-loop amplitudes for on-shell W~^ W~ production in the small-mass limit [72]. 
Starting from our results for the squared amplitude for W~^ W~ production (see 
Appendix A), we take the small-mass limit, namely rn^js —)■ 0 for fixed (t — m^')!s. 
Adjusting for overall conventions we found agreement with the results obtained in [72] 
in all contributions, except for F) {s,t) arising from the interference of two-loop 
diagrams in class C with the tree-level diagram in class A. From the discussion 
in [72], we could trace back this discrepancy to a different treatment of the vector- 
axial contributions in the fermionic loop in class C, resulting in a non-vanishing 
remainder even for zero-mass quarks. Since this appears to be inconsistent with 
charge parity conservation, we have good reasons to believe that the prescription 
used here as well as in [32] is the correct one (see our discussion in Section 4). 

6. Finally, we have compared numerically results both for the individual form factors 

Ej and for the full amplitudes and Mrll at tree-level, one-loop and two-loop 

order, with reference [32]. For the numerical evaluations of the helicity amplitudes 
we employed the package SOM [73]. We find full agreement with the results reported 
in [32], after a mistake in the calculation of one of the form factors was corrected in 
that reference. 

8 Numerical code and results 

We provide a C++ code for the numerical evaluation of the 9 finite form factors Ej for 
classes A, B and C. The implementation supports both, evaluation in the qT-scheme and 
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Figure 2. Real parts of the two-loop form factors for the process qq' V 1 V 2 in dependence 

of the relativistic velocity, /? 3 , and the cosine of the scattering angle, cos 6 * 3 , of the vector boson Vi. 
The virtualities of the vector bosons are set to = 2p|. 


in Catani’s original scheme. Further, it also provides the (alternative) 10 form factors Aj. 
The code is set up in form of a C++ library, which is supplemented by a simple command 
line interface. 

The code was optimised for speed and stability of the numerical evaluations, in par¬ 
ticular, by employing an appropriate functional basis for the multiple poly logarithms, see 
Section 5.2. We employ C++ templates to support evaluations with three different data 
types: double precision, quad precision and arbitrary precision using the CLN library [74]. 
The multiple polylogarithms are evaluated via their implementation [60] in the GiNaC li¬ 
brary [46], which also employs the CLN arbitrary precision capabilities. 

For the benchmark point of [32] no severe cancelations due to asymptotic kinematics 
take place. In this case our double precision implementation is accurate and gives at least 
10 significant digits for each of the Ej at the two-loop level. The evaluation of all Ej inch 
crossed variants, as needed for the physical amplitude, takes 150 ms on a single core of a 
standard desktop computer. Close to the phase space boundaries or in the high energy 
region, numerical cancelations lead to a significant loss of precision. In order to detect and 
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Figure 3. Real parts of the two-loop form factors Ej for the process qq' V 1 V 2 in dependence 
of the relativistic velocity, /? 3 , and the cosine of the scattering angle, cos 6 * 3 , of the vector boson Vi. 
The virtualities of the vector bosons are set to = 2p|. 

cure a possible instability, we compare the results obtained from evaluations with different 
precision settings and adaptively increase the precision until the target precision is met. We 
find the method to converge even in highly collinear configurations, where one needs to allow 
for a significant increase in the evaluation time though. Of course, also for unproblematic 
points in the bulk of the phase space, where the double precision results are actually 
accurate enough, our precision check requires additional run-time. For the aforementioned 
benchmark point we find an increase in the evaluation time to approximately 0.8 s on a 
single core. 

In Figures 2 and 3 we show numerical results for the class A and class C contributions 
to our 9 form factors Ej at the two-loop level. Note that these results were obtained with 
our C++ code and thus demonstrate the high numerical reliability of our implementation. 
We vary the relativistic velocity, = k/(s + p^ — p^}, and the cosine of the scattering 
angle, cos 6*3 = (2t -|- s — — p|)/ k, of the vector boson V. For the virtualities of the vector 

bosons we have set = 2 p|. All results are for Nf = 5 and given in the gT-scheme. The 
class A contributions in Figure 2 show pronounced structures in the collinear regions (see 
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(4.13) for the corresponding tree level coefficients). In contrast, the class C contributions 
in Figure 3 show no such features and are rather smooth functions in the full P^-cosO^ 
plane. 

9 Conclusions 

In this paper, we presented the derivation of the two-loop massless QCD corrections to the 
helicity amplitudes for massive vector boson pair production in quark-antiquark annihila¬ 
tion. The combination with leptonic decay currents allows to construct the two-loop QCD 
matrix elements relevant to four-lepton production. In this course, we computed all master 
integrals and optimised their representation for numerical performance. Our results ob¬ 
tained for the amplitudes provide a fully independent validation of a recent calculation [32]. 
We implemented our amplitudes in a C++ code for the fast and stable numerical evaluation 
of the amplitudes, which we provide together with our analytical results for public access at 
http://vvamp.hepforge.org. This opens up the path towards precision phenomenology 
in gauge boson pair production and improvements of the background predictions for Higgs 
boson studies and searches for physics beyond the Standard Model. 
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A Squared amplitudes for the on-shell production of vector-boson pairs 

In this Section we show how the general results described in this article can be used to obtain 
the squared amplitude for the process V 1 V 2 summed over spins and colours. Eor the 

calculations of the NNLO QGD corrections to on-shell ZZ [27] and W^W~ production [28] 
production, we directly computed the squared amplitudes using a dedicated setup based 
on our solutions for the equal-mass master integrals [25, 26]. We compared the results 
obtained in the two approaches and find full agreement. 

We denote the squared amplitude as 

{M\M) = T{s,t,pl,pl) = ^ |5';,i.(pi,p2,P3)e3(7’3)e4(P4)|\ (A.I) 

pol,colour 

which of course can be perturbatively expanded in powers of ag as 
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T{s,t,pl,pl) = (47ra) 


2 T^°\s,t,pl,pl) + T^^Hs,t,pl,pl) 


+ (^) T^‘^Hs,t,pl,pl) + Oial) 


(A.2) 


where we have 


= (A.3) 

T^^\s,t,plpl) = 2?li[{M^^^\M^^^)) , (A.4) 

T^‘^\s,t,plpl) = 23f? + (AlW|A^W). (A.5) 

It is easy to write a general expression of in terms of the coefficients 

and Aj™'^ or, equivalently, in terms of the Tj^'^ and Tj^\ simply by contracting the general 
decomposition (2.11) with itself and summing over colours and external polarisations us¬ 
ing (2.9). The result is quite involved and not particnlarly illuminating and we decided not 
to include it here explicitly. This general formula, in fact, is needed explicitly only in order 
to derive the 1-loopxl-loop corrections which can however also be easily ex¬ 

tracted from automated codes, and therefore we will not consider them here. On the other 
hand, if we limit ourselves to considering the contraction of the generic n-loop amplitude 
with the tree-level, i.e. m = 0, the results are much more compact. In the following two 
sections we will discuss the two explicit cases of on-shell ZZ and WW production, which 
were used for the calculations in [27, 28]. 


A.l The two-loop corrections to ZZ production 

In the case of gg —)• ZZ the tree-level is given by the two diagrams belonging to classes 
C = A, B. As far as two-loop corrections are concerned, the classes of diagrams that can 
contribute to ZZ production are C = A, B,C, see Section 4. By contracting the tree-level 
diagrams with the general amplitude (2.11) one easily finds 


N 






' {ZZ,{n)) (ZZ,{n)) 


nJZZ,{n)) {ZZ,{n))' 
Z Tj T;^q 


u t J 

(A.6) 

Igq are defined in (3.16). Each of the 
^(zz,(n)) obtained summing over the relevant classes of diagrams, re-weighted by 

appropriate coupling factors 


where N is the number of colours, while and R, 


r™) = + rf l’(") + TVzz rf l’(^) 


(A.7) 


[Cl 

where the rj ^ components of the tj are defined by a decomposition completely analogous 


to that for the Aj in (3.22), 


[(L|_)2 + (i?f )2] 

Nyy = ^ Nyy 

[(L|,)4 + 


(A.8) 
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and Nzz is defined in (4.1) . We have verified explicitly that as far at the tree-level and 
one-loop corrections are concerned, we have full agreement with the results in [69]. Similar 
but much more lengthy formulas can be derived for and we do not report 

them here for brevity. 


A.2 The two-loop corrections to W~^W production 

Let us consider now the case of qiQi —>■ W^W~, where the index i labels the flavour of the 
initial state quarks, qi = {u,d). At the tree-level, this process receives contributions from 
three diagrams, one in class A and the other two in class Fy, with V = Z, 7 . Let us start 
from the tree-level and one-loop corrections, where only diagrams in classes C = A, Fy can 
contribute. Following the notation of [70], we separate the contributions to the squared 
amplitude into three different form factors 




cfFf\s,t) 


„ts t{0) 




(0) 


(s,t) 


(A.9) 


23^ 



ww 


N 


c“ f/'> (s. t) - jfL. t)+c L'L.«) 


(A.IO) 


contains the squared contribution of diagrams in class C = A (i.e. diagrams where 
the production of the W^W~ pair is not mediated through a 7 or a Z boson). 
encapsulates instead the interference of the Fy -type diagrams (i.e. those where the W^W~ 
pair is produced via a 7 or a Z virtual boson) with diagrams in class C = A. Finally 
is given by the interference of the Fy-type diagrams with themselves. Again, following 
closely [70] we define then 


* 16 sin^ Ow ’ 

SS _ 1 ^ZW+W-{L^iq,+ Rqiq,) S \ f Czw+W-i^lq^ - R^q^) S 

' 2 s-rnlj 2 s-i 


(A.ll) 


where, as always, Cq. is the quark charge in units of e, with e > 0 , and the electroweak 
couplings , Rq^qi and Czw+W- are defined in (3.16) and (4.3). 

At two loops the decomposition (A.IO) must be enlarged since also diagrams belonging 
to class C = C start contributing to the amplitude. We therefore write the two-loop 
contribution as follows 


23? 


2 )' 


i,WW 


= N^cfFP{sA) + cP’**Ff'’(^^(s,t) 

-cf jf (s,t)-cP’*Vf]’(2)(s,t) + cfAf)(s,t)l , (A.12) 
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where we introduced the new couplings 


JC],« _ 
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J^C],ts _ 


32 sin^ 
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-N, 


9 ’ 


<^ZW+W- {^qiQi + 




4 s sin^ 


iV, 


s — nir. 


9 ■ 


(A.13) 


\C\ (2) [Cl (2) 

Here, the new form factors F- {s,t) and is,t) contain the contribution from the 

two-loop diagrams in class C = C. In deriving (A.13) we used the fact that for a fermion 
loop with an attached VH-pair we have 

^ww = \Y 1 = 4sin2g^ ^g ’ (A-14) 

qq' 


where Ng = Nf/2 is the number of generations of massless quarks running in the loop. 
Note that because of the flavour-change induced by the bosons, we limit ourselves to 
consider at most Nf = 4 massless quarks {u,d,c,s), i.e. two generations Ng = 2. Finally, 
the form factor K^\s,t) receives contributions only from one class of diagrams, C = Fy. 
At tree level we find that the different form factors can be obtained from 




2 r, 


[^],( 0 ) 


10 


— 4r- 
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f\s,t) =a(t. 


t 

[A,(o) , ^[A],(o) 
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+ T-r 


o ( ,[A,( 0 ) , ,[A],( 0 ) 
^ \ I a I /1 


[nm , .^[i"],(o) 


Af)(s,t) = 2(4 
At one loop and two loops we find instead 


- T, 


9 

imo) 


'10 
■r Ho 


'n A [A],(n)' 

^ Ho ~ ^ H 


2 ( 


Fl'^\s,t) = 23f? 

= 23f? 

K\^\s,t) = 23f? (Kf\s,t)F^^\s)) , 
and the two new form factors read 


(A.15) 

(A.16) 

(A.17) 


- T. 




(A.18) 


FP’^'^\s,t) = 23? ( - 


— 4t^ 


p],(2)' 


= 251? [2 - f 


(A.19) 

(A.20) 


where are the n-loop QCD corrections to the quark form factor. We have verihed 

that the tree-level and one-loop corrections, in the limit of equal virtualities of the massive 
vector bosons, agree with [70]. 
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B Schouten identities for the amplitude 


In this Appendix, we show how to reduce the number of independent form factors entering 
our helicity amplitudes by exploiting the 4-dimensionality of external states via Schouten 
identities. We document here a general way to derive such Schouten identities for the 
^RLL case. The LLL case proceeds in exactly the same way. 

We start off by fixing the helicities for a right-handed incoming quark current in the 
spinor helicity notation and we get 

SR{pi,P 2 ,P 3 ) = [ 2 1 ^ 3 !) + A2P^P2 +AsPip'^ + A4 pi^p2) 

+ [ 27 ^" 1 ) {A5P1 + AQP2) + [ 27 ^^ 1 ) {Ajp'l + Aspi^) 

+ Ag [2 1) + ^10 [2 1). 


As a hrst step we notice that we can collect [ 2^3 1 ) as an overall factor: 


[2^3 1)[M3 2) = TV 


1 + 75 

2 


j. 22 

= tu - P 3 P 4 . 


Multiplying and dividing by this allows to write the partonic amplitude as 
SR{pi,pt,P3) = [2j^3l) I + A 2 P 1 P 2 + Asp'^pi^ + Ai 


+ 

+ 


{A 5 P 1 + AQP 2 ) + 


32 ) [2 7 - 1 ) 


tu-plpl 


[ 1^3 2 ) [ 27 ^ 1 ) 

t u - pIpI 

, 2 127^3^ 1 ) + , 

tu-pipz tu-pipz 


{A 7 P 1 +ASP 2 ) 


[1^3 2) [ 27 ^ 7 " 1)}, (B.l) 


such that every spinor structure is a trace. We can then perform the traces recalling that 
the transversality of the leptonic decay currents allows to discard contributions proportional 
to P 3 or P 4 . In this way we get 


[1j^3 2) [ 27 ^ 1 ) - {u - pDpi - {t - pI)p'^ (B.2) 


and 


[1^3 2 ) [ 27 ^ 7 " 1 ) 


2 {u-pD + 2pi -{tu- plpl)g^^ 

2up'^P2 + 2pl p’ip'^ -2 {u-pI) p^pi 


(B.3) 


[1 ^3 2) [2 1 ) = -2{u- pI) _ 2^2 + 4 ^pi,p 3 ,P 2 ,p ^ 

-{tu- plpl)g^'' - 2 tpip^ + 2p|p^p^ - 2 (t-pl) , (B.4) 

where we introduced the Levi-Civita e tensor, with the following notation 


Moreover, note that the asymmetry between (B.3) and (B.4) is due to the transversality 
condition which effectively replaces P 3 —>• 0 and P 3 —>• Pi + P 2 • 
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Using (B.2),(B.3) and (B.4) we see that all 10 spinor structures can be written in terms 
of the following 11 structures: 




g 


f:Pi,P3,P2,M 

e Pi ; 


U u U u U u u 

PIPI, P 1 P 2 > P 2 P 1 > P 2 P 2 1 


e P 2 


fPl,P3,P2,l' 

e Pi . 

cPl’P2,^J■,^' 


fPl,P3,P2,l2 pP 


This does not appear to be any improvement with respect to the 10 structured we had 
before. It is nevertheless very easy to show that 2 out of these 11 structures can indeed be 
expressed as linear combinations of the remaining 9 by means of an anti-symmetrisation 
of the eP^P^ tensors. 

In order to see how this works in practice, we start off by considering gPi>P3,M.!^p2 • pi. 
By anti-symmetrising in 4 dimensions one easily finds 

^Pl,P3,P,<2 p^ . p^ ^ _^P3,P,y,P2 pi.pi- ^P,^,P2,P1 p^.pi- e^.P2,Pl,P3 pP _ ^P2,P1,P3,P pU 

which implies that eliminated by 


cPl^P3,P,’2 — 


2 


Pi-t 


^Pl,P2,P,l' _j_ ^Pl,P3,Pl,12 pP _ ^Pl,P3,P2,P pt^ 


(B.6) 


s V 2 

leaving us again with 10 structures. One more anti-symmetrisation can be used, namely 
consider eP3’P2,p,’2 p^ . where the momentum is defined as 

(u-pI\ ft- pI\ ^ ^ 

j + \~y~ j ^2 

such that r ■ Pi = 0, r ■ p 2 = 0- Proceeding as before we find 

fPl,P3,P2,P f. 

fc r, 9n ,, . Hi t 


f:Pl,P2,P,12 _ 


t u - pIpI 


[(tt — p 1) P2 + {t — p\) Pi] -i- 


pi,P3,P2,l^ 


tu-pjpl 


[{u - pI) Pi + {t - pI) pf^] . 

(B.7) 


It becomes clear that using these two relations we can eliminate completely ePi’P 2 .M.i^ 
and ePi’P 3 ,M ,!2 favour of the remaining 9 structures. In particular these relations can be 
rephrased in terms of the original spinors in (B.l) giving two Schouten identities for the 
spinor lines: 


[1 ^3 2) [2 Yhl'' 1) = {tu- pIpI) 


-{ptp^2-PiP^2)-9^'' 


1 


+ -[(“- pi) Pi - (i - pi) P2] [1 h 2) [2 Y 1) 

- ^ [(n - s - pi) p^ + (n - pi) p^j [^3 2 ) [2 1 ), 


(B.8) 


[1 ^3 2) [2 1) = {tu- P3P4) 


{P 1 P 2 - P 1 P 2 ) - g'"‘' 


1 


+ - [{t-pl)P 2 - {u-pI)p^] [1^32) [27^1) 

- ^ [(i - p 1) K + (t - ^ - pi) P2] [1^3 2) [2 1). 


(B.9) 
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The corresponding relations for the spinors of the left-handed partonic currents can 
be found by simply permuting pi -H- p 2 . Using (B.8),(B.9), and the corresponding ones 
for the left-handed partonic current, we eliminate 2 of the structures in (B.l) in favour of 
g^'^, plus the remaining 8 structures, and then proceed by contracting with the left-handed 
leptonic decay currents (3.8). As a result one easily arrives at formulae (3.19) and (3.20). 


C Conversion to Catani’s original IR snbtraction scheme 


In Section 6 we derived the finite remainder of the one- and two-loop helicity amplitude 
coefficients 12 in a subtraction scheme which is particularly well-suited for qx subtrac¬ 
tion [68]. In this Appendix we show how these results can be converted to Catani’s original 
scheme [67]. Starting from the UV-renormalised coefficients defined in (6.5) at renormali¬ 
sation scale p?, we write the finite remainders in Catani’s scheme as 

= n'"’ - 'f {') n'” - ‘S(‘) . (O-i) 


where Catani’s subtraction operators are defined as 

0^7 /i q\/ ,.2\^ 


= -Cf 


1 3 

r(l-e) V? 
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^2 (f) — (f) + ) + 


rCi 


2 / 3 o 


e-^7r(l - 2e) fl3o 


r(l-e) 


— +k) If(2e) + Il(2)(g) (0 2) 


with 




and since a qq pair is the only coloured state we have 

pE 7 / ,. 2^26 


(C.3) 


Il( 2 )(g) = 


4er(l - e) 
.2 
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X 2Cf 


IT 


13 


245 23 


— - 6 Ca - - 1 Cf + ( —Cs + — - —vr^ j Ca + ( — - — 

n sd o ; ^ -r i o'’■i ^ 216 48 / V 12 '''■ 


25 


TpNf . 

(C.4) 


In this article, we present our results for = s. Note that upon expansion in e both if (e) 
and if (e) generate imaginary parts whose sign is fixed by the prescription s —?• s -|- i 0^ . 

By comparing (6.6) with (C.l) one can show that the parts of the finite, complex form 
factors in Catani’s original scheme [67], can be obtained from those in the gF-scheme [68] 
according to 


q( 1),finite _ q( 1),finite i A r (->(0),finite 

^ ^Catani ~ ^ ^ ^ ’ 

(2),finite _ q( 2),finite , A r q( 1),finite i A r q( 0),finite 
“Catani ~ “iJt , 


(C.5) 
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with the finite scheme conversion coefficients given by 


All = Cf 

' 162 


(C.6) 


187 


/961 11 


13 


432 


vr + v;7rC3 + —vr - —Cs 


9 2 , 1 4 


72 
3 5 


96 


V216 72 


+ +- 7 r +^’^(^3-4^ +6C3 

.r ^ ('41 97 2 17, . / 65 1 2 

1 ' 81 216 36^ V 108 36 


(C.7) 


where we have set = s to match the convention for our final results. Notice that, in order 
to obtain the finite remainders of the two-loop amplitudes in the two different schemes, only 
the finite pieces of the latter are required, and in particular the 0{e) terms of the one-loop 
amplitudes are not needed, as expected. Note, moreover, that the conversion coefficients 
are complex, due to the fact that the original formulation of IR subtraction [67] factored 
out a phase for time-like pairs of partons from both the collinear and soft contributions, 
while in the gr-scheme [68] this phase factor is associated only with the soft contributions, 
in line with the structure of IR factorisation [77, 78] at higher loop order. 
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